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INTRODUCTION 

Ever  since  their  discovery  (e.g.,  Brown  and  Roshko,  1974),  coherent  structures  in 
turbulent  shear  flows  have  proven  to  be  important  in  enhancing  the  mixing  and  entrainment  of 
chemical  species  into  the  turbulent  region.  An  understanding  of  the  evolution  of  coherent 
structures  is  necessary  to  predict  the  mixing  rates  in  turbulent  jets  accurately.  Furthermore, 
their  importance  in  reacting  flows  becomes  apparent  if  the  large-scale  structures  could  be 
enhanced  to  increase  mixing.  In  the  past,  both  active  and  passive  techniques  aimed  at 
controlling  the  evolution  of  coherent  structures  have  been  used.  By  enhancing  the  spreading 
rate  in  the  shear  layer,  high  mixing  rates  between  the  fuel  and  the  air  can  be  achieved.  This  is 
owing  to  the  simple  observation  that  entrainment  into  the  jet  cannot  occur  faster  than  the  jet 
growth.  Obviously,  increased  mixing  between  the  fuel  and  the  air  streams  would  also  result  in 
enhanced  combustion. 

Of  the  active  methods,  acoustic  forcing  of  the  shear  layer  at  various  frequencies  has  been 
shown  to  enhance  the  generation  of  highly  coherent  large-scale  structures  (e.g..  Ho  and  Huang, 
1982;  Oster  and  Wygnanski,  1982).  Studies  have  shown  that  the  spreading  rate  can  be 
dramatically  increased  when  acoustic  forcing  is  applied  at  the  subharmonic  of  the  most  unstable 
frequency.  This  can  result  in  the  simultaneous  merging  of  several  vortices,  even  when  the 
forcing  amplitude  is  small.  Experiments  by  Winant  and  Browand  (1974)  have  shown  that  the 
vortex  rollup  and  pairing  process  is  the  dominant  mechanism  of  shear  layer  growth.  Com¬ 
bustion  enhancement  in  coaxial  jets  of  fuel  and  air  was  demonstrated  by  Gutmark  et  al.  (1989) 
by  forcing  both  the  jets  at  prescribed  frequencies  that  were  related  to  the  most  unstable  modes. 

Passive  control  techniques  have  also  been  extensively  studied.  Ho  and  Gutmark  (1987), 
Gutmark  and  Ho  (1985),  Hussain  and  Husain  (1989)  and  Schadow  et  al.  (1986)  have  ex¬ 
perimentally  studied  the  entrainment  characteristics  of  jets  with  different  cross  sections.  It  was 
observed  (Ho  and  Gutmark,  1987)  that  an  elliptical  jet  with  an  aspect  ratio  of  2:1  was  about  3 
to  8  times  more  effective  in  entraining  fluid  than  a  circular  jet  with  the  same  cross-sectional 
area.  The  reason  for  this  mixing  enhancement  was  attributed  to  a  complex  series  of  distortions 
of  the  coherent  elliptical  vortex  rings  due  to  their  self-induced  flow  and  to  their  mutual  inter¬ 
actions.  These  nonlinear  interactions  result  in  the  phenomenon  of  "axis  switching,"  where  the 
orientation  of  the  major  and  minor  axes  of  the  elliptical  vortex  ring  appears  to  switch  during  its 
propagation.  The  number  of  switchings  and  the  locations  are  strongly  dependent  upon  initial 
conditions,  the  aspect  ratio  and,  when  excited,  the  Strouhal  number  and  the  excitation  level 
(Hussain  and  Husain,  1989).  The  deformation  and  axis  switching  of  the  elliptical  jets  appear  to 
play  a  major  role  in  the  increased  entrainment  observed  with  elliptical  jets  (Hussain  and  Husain, 
1989).  How  this  entrainment  process  is  modified  when  there  is  heat  release  is  not  clearly  known 
at  present.  In  a  chemically  reacting  flow  when  heat  is  released  the  gas  density  is  lower  in  the 
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reacting  zone.  The  dynamics  of  these  elliptical  structures  can  be  affected  by  this  density 
variation. 

To  understand  the  complex  flow  features  associated  with  the  propagation  and  breakdown 
of  elliptical  jets  with  and  without  heat  release,  a  numerical  study  was  undertaken.  Due  to  the 
change  in  affiliation  of  the  principal  investigator.  Dr.  Wen-Huei  Jou,  this  research  was 
terminated  after  6  months  of  study.  This  report  documents  the  findings  of  the  research  in  that 
period  and  presents  some  preliminary  results  of  a  study  that  was  carried  out  subsequent  to  the 
termination  of  the  contract. 


TECHNICAL  OBJECTIVES  OF  THIS  STUDY 

A  research  effort  involving  numerical  simulations  was  undertaken  in  order  to  understand 
the  details  of  the  mixing  processes  associated  with  elliptical  jets  with  and  without  heat  release. 
The  effort  was  intended  to  study  the  details,  using  full  three-dimensional  simulations,  of  the 
development  of  the  nonlinear  instability  waves  (Koshigoe  et  al.,  1988),  the  interaction  (merging) 
of  elliptical  vortex  rings  and  the  associated  entrainment  process.  The  effect  of  heat  release  on 
the  evolution  of  the  mixing  region  was  also  to  be  investigated.  Flow  field  statistics  were  to  be 
computed  throughout  the  development  of  the  jet  to  characterize  the  mixing  properties  of  the 
simulated  flow  field. 


SIGNIFICANT  RESULTS  IN  THE  PAST  YEAR 

In  the  past  year,  the  research  effort  was  entirely  directed  towards  constructing  a  fourth- 
order  accurate,  compact  finite-difference  scheme  for  chemically  reacting  flows.  Earlier  studies 
of  chemically  reacting  flows  have  employed  numerical  schemes  using  pseudospectral  techniques 
(McMurtry  et  al.,  1986,  1989).  These  simulations  were  restricted  to  spatially  periodic  flows. 
The  study  of  more  realistic  flows,  such  as  in  an  elliptical  jet,  requires  a  numerical  simulation  of 
spatially  evolving  flows,  due  to  the  decomposition  of  the  initial  disturbance  into  a  set  of 
complex  modes  with  varying  dispersive  properties.  In  addition,  to  resolve  the  shear  layer 
regions,  a  stretched  nonuniform  grid  in  the  cross-sectional  plane  is  required  to  prevent 
boundary  effects  in  simulating  axis  switching. 

For  these  reasons,  a  compact  higher-order  scheme  is  more  flexible  in  treating  the  bound¬ 
ary  conditions.  There  are  many  advantages  of  a  compact  scheme.  First,  a  higher  order  of 
accuracy  (fourth  order  in  this  case)  can  be  achieved  using  a  compact  stencil  of  only  three  grid 
points.  This  high  accuracy  is  obtained  by  solving  for  the  derivatives  simultaneously  at  three 
points  using  a  tridiagonal  solver.  The  compact  scheme  also  has  a  smallei  truncation  error  and 
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smaller  phase  error  than  conventional  centered,  fourth-order  finite-difference  approximations. 
These  properties  are  very  important  when  performing  simulations  of  the  spatial  development  of 
unsteady  flows.  In  the  following,  we  briefly  describe  a  simple  test  problem  that  was  used  to 
study  the  effectiveness  of  the  compact  scheme  to  resolve  wave  motion  without  significant  phase 
error. 


COMPACT  SCHEMES 

Consider  a  one-dimensional  unsteady  wave  equation  of  the  form 

du/dt  =  du/dx  (i) 

which  is  to  be  solved  subject  to  some  specified  initial  condition  at  /  =  0.  A  fourth-order 
compact  scheme  approximation  of  the  first  differential  spatial  operator  can  be  written  down  as 

du/dx  =  [1  +  62/ 6]1  [(u/+,  -  ui-t  )/2Ax]  (2) 

where  b2  is  the  central  second  difference  spatial  operator,  and  the  discretization  is  on  a  uniform 
grid.  Extension  to  a  nonuniform  grid  is  straightforward.  This  spatial  discretization  can  be 
combined  with  time-accurate  schemes,  such  as  the  second-order  Adams-Bashforth  scheme  or 
the  Runge-Kutta  scheme,  to  obtain  a  high-order  spatially  and  temporally  accurate  numerical 
algorithm. 

A  detailed  study  using  the  one-dimensional  wave  equation  (1)  subject  to  an  initial 
disturbance  of  the  form  sin  ( kx )  at  /  =  0  was  carried  out  to  evaluate  the  various  compact 
schemes.  The  length  of  the  one-dimensional  duct  L  was  chosen  to  be  unity,  and  256  grid  points 
were  used  to  discretize  the  flow  field.  Periodic  boundary  conditions  were  used  at  x  =  0  and 
x  =  L.  The  nondimensional  wavenumber  k  of  the  initial  disturbance  was  varied  to  study  the 
effect  of  resolution  of  the  wave.  Note  that  an  increase  in  the  wavenumber  implies  a  decrease  in 
the  wavelength  of  the  disturbance.  Thus,  for  a  fixed  number  of  grid  points,  an  increase  in  k 
results  in  a  waveform  that  is  resolved  using  a  lesser  number  of  grid  points.  For  example,  for  a 
wavenumber  of  k  =  40,  the  wave  is  resolved  using  around  six  grid  points. 

The  primary  objective  of  this  numerical  experiment  was  to  determine  the  error  in  the 
phase  speed  of  the  initial  disturbance  during  the  temporal  evolution  of  the  wave.  The  fourth- 
order  and  sixth-order  compact  schemes  were  evaluated  using  time-marching  methods  such  as 
Adams-Bashforth  and  Runge-Kutta.  The  higher-order  accuracy  of  the  compact  scheme  in 
predicting  the  phase  speed  of  a  simple  wave  was  demonstrated  in  this  study.  Some  pertinent 
results  are  shown  here. 
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Figure  1  shows  the  error  in  the  phase  as  a  function  of  wavenumber  for  various  time- 
marching  schemes  for  the  fourth-order  compact  scheme  with  third-order  boundary  conditions. 
It  can  be  seen  that  the  Adams-Bashforth  scheme  with  a  CFL  of  0.1  can  predict  the  phase  speed 
of  the  wave  within  0.1  percent  for  the  high  wavenumber  k  =  40.  Although  Runge-Kutta 
schemes  could  be  used  with  a  higher  CFL,  the  overall  error  in  the  phase  speed  exceeds  0.1 
percent  by  k  »  30.  The  wave  motion  in  a  portion  of  the  duct  can  be  seen  in  Figure  2,  which 
shows  the  shape  of  the  wave  for  the  wavenumber  k  =  40  at  three  different  times.  It  can  be 
seen  that  the  wave  is  adequately  resolved  using  only  six  grid  points. 

Figure  3  shows  a  comparison  of  the  fourth-order  [0(4)]  and  the  sixth-order  [0(6)] 
compact  schemes  using  the  Adams-Bashforth  time-stepping  scheme.  It  appears  that,  for  the 
same  wavenumber,  the  fourth-order  compact  scheme  is  superior  to  the  sixth-order  scheme. 
Based  on  this  study,  the  decision  was  made  to  implement  a  fourth-order  compact  scheme  in  the 
three-dimensional  solver. 


THREE-DIMENSIONAL  NAYTER -STOKES  CODE  DEVELOPMENT 

A  new  hybrid  code  was  developed  in  this  study  that  uses  a  fourth-order  compact  scheme 
in  the  streamwise  direction  and  pseudospectral  methods  in  the  other  two  directions.  In  spectral 
methods,  the  solution  of  a  differential  equation  is  represented  as  a  truncated  series  of  known 
orthogonal  functions.  If  the  solution  to  a  problem  is  infinitely  differentiable,  the  error  in  a 
properly  implemented  spectral  method  can  be  shown  to  go  to  zero  faster  than  any  finite  power 
of  N  as  N  gets  larger,  where  N  is  the  number  of  modes  retained  in  the  series  expansion.  Other 
advantages  of  spectral  methods  over  more  commonly  used  finite-difference  techniques  are  that 
the  numerical  diffusion  and  dispersion  (phase)  cirors  are  very  small.  In  applications  to  reacting 
flows,  where  steep  gradients  of  the  chemical  reactants  can  result  and  where  reaction  rates  are 
diffusion  controlled,  numerical  diffusion  and  dispersion  can  be  major  sources  of  error. 

The  main  shortcoming  of  the  spectral  method  is  due  to  its  limited  application  in  problems 
with  complex  geometries.  The  choice  of  the  expansion  function  is  determined  by  two  factors: 
(1)  the  physical  boundary  conditions  of  the  problem  and  (2)  the  existence  of  fast  transform 
techniques  to  sum  the  truncated  series  efficiently.  In  a  spectral  method,  the  boundary  condi¬ 
tions  of  the  series  expansion  must  be  consistent  with  the  physical  boundary  conditions.  In  the 
application  of  Fourier  series  (the  most  common  expansion  functions  used  in  spectral  methods), 
we  are  limited  to  studying  problems  with  periodic,  even  or  odd  symmetries  at  the  boundaries. 

In  this  study  we  are  interested  in  the  spatial  development  of  the  elliptical  jet.  Thus,  the 
inflow  and  outflow  boundary  conditions  associated  with  streamwise  development  must  be 
implemented.  This  is  very  difficult  using  spectral  methods  but  can  be  accomplished  relatively 
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easily  using  the  compact  scheme  described  above.  As  shown  above,  the  phase  errors  associated 
with  the  fourth-order  compact  scheme  are  quite  low,  and  extension  of  the  method  to  a 
nonuniform  grid  is  straightforward. 

A  three-dimensional  Navier-Stokes  solver  using  spectral  methods  in  the  transverse  and 
spanwise  directions  and  the  fourth-order  compact  scheme  in  the  streamwise  direction  has  been 
developed.  All  the  subroutines  for  this  code  are  written  in  FORTRAN,  allowing  it  to  be 
portable  to  any  machine  with  minimal  modifications.  The  development  of  the  basic  code  was 
completed,  and  the  debugging  of  most  of  its  components  has  also  been  completed. 


RECENT  RESULTS 

Subsequent  to  the  termination  of  this  study,  the  development  of  the  numerical  code  was 
continued  by  Dr.  Patrick  McMurtry,  who  was  involved  in  the  earlier  development  work.  Some 
recent  results  of  the  initial  tests  are  described  here. 

The  basic  issue  in  implementing  the  compact  scheme  in  a  three-dimensional  solver  is  the 
proper  implementation  of  the  inflow  and  outflow  boundary  conditions.  To  study  this  feature, 
the  three-dimensional  code  was  reduced  to  a  two-dimensional  solver  that  used  spectral 
techniques  in  the  normal  direction  and  the  fourth-order  compact  scheme  in  the  streamwise 
direction.  A  two-dimensional  planar  jet  has  been  simulated  using  this  technique. 

The  flow  field  was  initialized  with  a  free-stream  velocity  in  the  top  half  of  the  domain 
(Uj)  with  a  value  of  1  and  in  the  bottom  half  ( U2 )  with  a  value  of  zero.  The  Reynolds  number 
based  on  the  velocity  difference  across  the  layer  was  200,  and  the  grid  resolution  was  100  x  32 
(x,  y).  At  the  inlet,  a  periodic  perturbation  was  added.  This  perturbation  corresponds  to  the 
most  unstable  wavelength  of  a  hyperbolic  tangent  mean  velocity  profile.  At  the  outflow,  a 
convective  condition  was  applied.  Mathematically,  this  condition  can  be  expressed  as 

d<t>/dl  +  Uc  df/dx  =  0  (3) 

where  Uc  is  a  convective  velocity  and  is  taken  to  be  (C/ 7  +  U2)/ 2. 

Figure  4  shows  the  evolution  of  the  spatially  developing  shear  layer  for  the  case  when 
only  the  fundamental  perturbation  was  used  to  initialize  the  flow.  As  expected,  after  the  shear 
layer  rolls  up  and  the  fundamental  reaches  a  maximum  amplitude,  the  vortex  is  convected  out 
without  any  further  growth. 

Figure  5  shows  the  case  in  which  the  subharmonic  perturbation  was  also  added.  The 
development  of  the  unstable  modes  is  clearly  evident.  A  pairing  process  between  the  two 
neighboring  vortices  is  beginning  near  the  outflow.  A  much  larger  computational  domain  will 
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be  needed  to  simulate  the  complete  pairing  process.  Note  that  at  the  outflow,  the  structures  are 
convected  out  without  affecting  the  internal  flow  field. 

These  simulations  clearly  demonstrate  that  the  basic  issues  involved  in  implementing  the 
compact  scheme  in  a  full  Navier-Stokes  solver  have  been  successfully  resolved.  Making  the 
current  two-dimensional  solver  into  a  full  three-dimensional  solver  should  not  pose  any  major 
problems,  since  all  the  elements  of  the  three-dimensional  solver  are  already  present  and  were 
merely  being  suppressed  to  simulate  the  two-dimensional  flow. 


RECOMMENDATIONS  FOR  FUTURE  EFFORT 

At  this  stage,  the  three-dimensional  code  development  is  mostly  completed  with  the 
resolution  of  the  issues  related  to  the  outflow  and  inflow  boundary  conditions.  With  some 
minor  effort,  this  code  can  be  made  ready  to  carry  out  a  serious  investigation  of  the  instability 
of  jets  with  and  without  external  excitation.  In  particular,  the  evolution  of  an  elliptical  jet 
subject  to  initial  disturbance  can  be  simulated.  Subsequently,  the  effects  of  heat  release  can  be 
included  v  evaluate  how  ’he  mixing  and  entrainment  processes  enhance  comb  union  efficiency. 
Comparisons  with  experimental  data  and  other  numerical  studies  could  be  made  to  validate  and 
demonstrate  the  capability  of  this  code  which  is  now  in  a  near-operational  status. 
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Wavenumber,  k 


Figure  1.  Phase  error  during  wave  motion  in  a  one-dimensional  duct  as  a  function  of 
wavenumber  using  different  time-stepping  schemes  and  the  fourth-order 
compact  scheme. 


A*B:  Adams-Bashforth  scheme  0(A i2) 
R*K:  Runge-Kutta  scheme  0(A i4) 
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Figure  2.  The  waveform  in  the  duct  at  three  different  times.  The  wa\e  is  rcsohed  b>  si\ 
grid  points  using  the  fourth-order  compact  scheme. 


TR-494/12-89 


9 


%  Error 


2  4  6  8  10  12  14  16  18  20  22  24  26  28  30 

Wavenumber,  k 


Figure  3. 


Phase  error  during  wau  motion  in  a  one-dimensional  duct  as  a  function  of 
wa\enumber.  Comparison  of  the  fourth-order  and  the  sixth-order  compact 
schemes. 
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Figure  5.  Vortex  rollup  and  pairing  in  a  spatially  developing  two-dimensional  shear  la>er 
using  the  hybrid  scheme.  Shear  layer  perturbed  by  the  fundamental  and 
subharmonic  modes. 
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